Traditionally, FEM or Finite Element Analysis (FEA) is one of the most popular computer aided engineering tool for engineers and scientists to model and solve engineering problems relating to complex systems such as three-dimensional non-linear structural design and analysis. FEA derives its name from the manner in which the geometry of the object under consideration is specified. With the advent of the modern digital computer, FEA has been implemented as FEA software. Basically, the FEA software is provided with a model of the geometric description and the associated material properties at each point within the model. In this model, the geometry of the system under analysis is represented by solids, shells and beams of various sizes, which are called elements. The vertices of the elements are referred to as nodes. The model is comprised of a finite number of elements, which are assigned a material name to associate with material properties. The model thus represents the physical space occupied by the object under analysis along with its immediate surroundings. The FEA software then refers to a table in which the properties (e.g., stress-strain constitutive equation, Young's modulus, Poisson's ratio, thermo-conductivity) of each material type are tabulated. Additionally, the conditions at the boundary of the object (i.e., loadings, physical constraints, etc.) are specified. In this fashion a model of the object and its environment is created.
However, FEM in certain circumstances has its drawbacks. For instance, FEM may not be advantageous due to numerical compatibility condition is not the same as the physical compatibility condition of a continuum. In particular, in Lagrangian type of computations, one may experience mesh distortion, which can either end the computation altogether or result in dramatic deterioration of accuracy. In addition, the FEM often requires a very fine mesh in problems with high gradients or a distinct local character, which can be computationally expensive. In a procedure called Arbitrary Lagrangian Eulerian (ALE) formulations, its objective is to move the mesh independently of the material so that the mesh distortion can be minimized. But the distortion in general still creates several numerical errors for large strain and/or high speed impact simulations. Furthermore, a mesh may carry inherent bias in numerical simulations, and its presence becomes a nuisance in computations. An example is the strain localization problem, which is notorious for its mesh alignment sensitivity. Therefore it is computationally efficacious to discretize a continuum by a set of nodal points without mesh constraints.
As a result, in recent years, meshfree method has also been used to solve these drawbacks suffered in FEM. However, meshfree method suffers its own problems. For example, the computational costs are generally too high comparing to that of FEM. It is very difficult and computational expensive to impose essential boundary conditions in meshfree method.
Neither FEM nor meshfree method can numerically predict every possible structural behaviors under various environmental loads, for example, a numerical phenomenon called volumetric locking for simulating structural behaviors in compressible and/or near-incomprssible regions (e.g., rubber, material under large plastic deformation). FIG. 1A depicts a simplified two-dimensional example illustrating the volumetric locking effect. Element 102 is fixed in two of its four sides. As a result, node 112 is fixed in space thus not moving in numerical simulation. However, this does not represent a physical reality.
Volumetric locking can be described using a set of mathematical derivation as follows:
Consider a plane-homogeneous isotropic linear material body which occupies a bounded polygonal domain Ω in 2 with boundary Γ=∂Ω. The problem is restricted to small deformation and the infinitesimal strain tensor g is defined as a function of the displacement u by
                              ɛ          ⁡                      (            u            )                          =                              1            2                    ⁢                      (                                          ∇                u                            +                              ∇                                  u                  T                                                      )                                              (        1        )            For a prescribed body force f, the governing equilibrium equation in Ω reads∇·σ+f=0  (2)where σ is the symmetric Cauchy stress tensor. With the fourth-order elasticity tensor denoted by C, the constitutive equation is given byσ=Cε=2με+λ(trε)I in Ω  (3)where I is the identity tensor of order four. Symbols μ and λ are Lamé constants which are related to the Young's modulus E and Poisson ratio v by
                              μ          =                      E                          2              ⁢                              (                                  1                  +                  v                                )                                                    ,                  λ          =                      vE                                          (                                  1                  +                  v                                )                            ⁢                              (                                  1                  -                                      2                    ⁢                    v                                                  )                                                                        (        4        )            
For simplicity, the displacement is assumed to satisfy the homogeneous Dirichlet boundary condition u=0 on ΓD and the traction t is distributed over the Neumann boundary ΓN with Γ=ΓD∪ΓN and ΓD∪ΓN=0. The variational form of this problem is to find the displacement uεV such thatA(u,v)=l(v)∀vεV  (5)where the space V=H01(Ω) consists of functions in Sobolev space H1(Ω) which vanish on the boundary and is defined byV(Ω)={v:vεH1,v=0 on ΓD}  (6)
The bilinear form A(.,.) and linear functional l(.) in Eq. (5) are defined by
                                                                        A                :                VxV                            ->                                                                          A                ⁢                                  (                                      u                    ,                    v                                    )                                            =                                                2                  ⁢                                                                          ⁢                  μ                  ⁢                                                            ∫                      Ω                                                                                                            ⁢                                          ɛ                      ⁡                                              (                        u                        )                                                                                            :                                                                            ɛ                      ⁡                                              (                        v                        )                                                              ⁢                                                                                  ⁢                                          ⅆ                      Ω                                                        +                                      λ                    ⁢                                                                  ∫                        Ω                                                                                                                      ⁢                                                                        (                                                      ∇                                                          ·                              u                                                                                )                                                ⁢                                                  (                                                      ∇                                                          ·                              v                                                                                )                                                ⁢                                                                                                  ⁢                                                  ⅆ                          Ω                                                                                                                                                                            (        7        )                                                                                    l                :                V                            ->              •                                                                          l                ⁡                                  (                  v                  )                                            =                                                                    ∫                    Ω                                                                                                  ⁢                                                            f                      ·                      v                                        ⁢                                                                                  ⁢                                          ⅆ                                                                                          ⁢                      Ω                                                                      +                                                      ∫                                          Γ                      N                                                                                                                      ⁢                                                            t                      ·                      v                                        ⁢                                                                                  ⁢                                          ⅆ                      Γ                                                                                                                              (        8        )            
The bilinear form A(.,.) in the above equation is symmetric, V-elliptic and continuous on V. By Korn's inequality and Lax-Milgram theorem, there exists a unique solution to the problem.
The standard Galerkin method is then formulated on a finite dimensional subspace Vh⊂V employing the variational formulation of Eq. (5) to find uhεVh such thatA(uh,vh)=l(vh)∀vhεVh  (9)
The second term on the right hand side (RHS) of Eq. (7) resembles the penalty term similar to the classical penalty method in the Stokes problems. For the energy in Eq. (7) to remain finite as λ→∞ (or v→0.5), the following constraint must be enforced∇·u=0 for uεV  (10)Similarly in the finite dimensional space, we have∇·uh=0 for uhεVh  (11)
The solution of Eq. (9) using a low-order approximation for displacement field in general does not satisfy the constraint in Eq. (11) and fails to converge uniformly as v→0.5. This is the volumetric locking in the near-incompressible problem. In other words, volumetric locking occurs when the approximation space Vh is not rich enough for the approximation to verify the divergence-free condition ∇·uh=0. Shown in FIG. 1B are eight non-physical modes for 3×3 discretization using EFG—a prior art meshfree procedure.
An eigenvalue analysis using the left hand side (LHS) of Eq. (9) can be employed to further reveal the locking phenomena in the low-order displacement-based finite element and meshfree methods. The eigenvalue analysis of one rectangular quadrilateral bilinear (Q1) finite element with 2×2 Gauss integration points is analyzed first and two non-physical locking modes are obtained and shown in FIGS. 1C and 1D. The eigenvalues grow unbounded as v approaches to 0.5, which is not expected physically. Accordingly, eight non-physical locking modes are obtained from the eigenvalue analysis using four rectangular Q1 elements. Similar to the results in the four Q1 finite element model, the same number of non-physical locking modes are generated by Element-Free Galerkin (EFG) method as displayed in FIG. 1B using the first-order approximation with normalized nodal support equal to 1.5. The increase of nodal support and order of Gaussian quadrature rule does not completely eliminate the volumetric locking in EFG method. The number of non-physical modes remains unchanged.
Based on the aforementioned drawbacks, shortcomings and problems, it would, therefore, be desirable to have an improved computer-implemented method for numerically simulating structural behaviors in compressible and/or near-incompressible regions.